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ABSTRACT 

Self similarity allows for analytic or semi-analytic solutions to many hydro- 
dynamics problems. Most of these solutions are one dimensional. Using linear 
perturbation theory, expanded around such a one-dimensional solution, we find 
self-similar hydrodynamic solutions that are two- or three-dimensional. Since 
the deviation from a one-dimensional solution is small, we call these slightly two- 
dimensional and slightly three-dimensional self-similar solutions, respectively. As 
an example, we treat strong spherical explosions of the second type. A strong 
explosion propagates into an ideal gas with negligible temperature and density 
profile of the form p(r,0,4>) = r~ w [l + aF(6,<j))], where u > 3 and ff < 1. An- 
alytical solutions are obtained by expanding the arbitrary function F(9, 0) in 
spherical harmonics. We compare our results with two dimensional numerical 
simulations, and find good agreement. 

Subject headings: hydrodynamics, shock waves, instabilities 
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Introduction 



Astrophysics supplies ample examples of hydro dynamic problems that admit self-similar 



solutions. In supernovae explosions (Koo and McKee 1990 Chevalier 1976) a shock wave 



is created by the release of an immense amount of energy during a short time in the center 
of an exploding star. When the shock wave propagates into the surrounding medium, the 



hydrodynamics is described by the the Sedov- Taylor solutions (Sedov 1946; Von-Neumann 



1947 Taylor 1950 Waxman and Shvarts 1993). Gamma-ray bursts provide a relativistic 



analog of that (Blandford and McKee 1976 Best and Sari 2000; Sari 2006 Pan and Sari 



2006). If the external medium is spherical, these are one-dimensional solutions. However, if 



the external density has angular dependence, it will cause the shape of the shock, and the 
flow behind it, to deviate from sphericity. 

An inherently two-dimensional version of this problem is the explosion in half space. 
Here, space is assumed to be empty on one side of a plane, while the other side is filled with 
an ideal gas with constant density. A large amount of energy is then released at a point on 
the surface. This describes the propagation of Shockwaves in the process of cratering caused 
by large impacts on a planetary surface. Qualitatively, this problem and its self-similar 



nature was described by Zel'Dovich and Raizer (1967), but a two-dimensional self-similar 



solution was not developed there. 

Here, we obtain two-dimensional and three-dimensional self-similar solutions that 
deviate only slightly from some known one-dimensional solution. We show that when 
treating such solutions as perturbations, the analysis is analogous to the treatment of 



stability (Ryu and Vishniac 1987 Goodman 1990 Chevalier 1990 Sari et al. 2000 Kushnir 



et al. 2005). We call these solutions slightly two-dimensional or slightly three-dimensional 



self-similar solutions. As a working example, we analyze small deviations from sphericity 
in the case of the strong explosion problem with external density falling as a power law of 
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distance p oc r w , where a; > 3. Solutions with these values of u are known to be self-similar 



solutions of the second type (Waxman and Shvarts 1993) 



In Sec. [2] we briefly review the main features of the one- dimensional solution 
which serves as the unperturbed solution for our analysis. In Sec. [3j we discuss the 
perturbation formalism for this problem and find the slightly two-dimensional and slightly 
three-dimensional self-similar solutions for density perturbations proportional to a spherical 
harmonic. These solutions are then demonstrated using the values of u and 7 for which the 
unperturbed solution is analytic in Sec. |4j The solution is then given for small deviations 
from sphericity with arbitrary angular dependence in Sec. |5j while in Sec. [6] the small I 
limit is investigated. Our semi-analytic solutions are then favorably compared with full 
fluid-dynamic simulations in Sec. [7| Finally, in Sec. [8] we give our concluding remarks. 



2. The One-Dimensional Self-Similar Solution 



Here we summarize the formalism leading to the one-dimensional self-similar solution 



(Waxman and Shvarts 1993). The discussion here follows Sari et al. (2000). Consider the 



Strong Explosion Problem in which a large amount of energy is released at the center of 
a sphere of ideal gas with a density profile decreasing with the distance from the origin 
according to p = Kr~ u , forming a strong outgoing shock wave. 

This problem was first investigated by Sedov (1946), Von-Neumann (1947), and Taylor 



(1950), who found the solutions for 00 < 5, known as the Sedov- Taylor solutions. Waxman 



and Shvarts (1993) showed that the Sedov- Taylor solutions are valid only for tu < 3, where 



the solutions are known as self-similar solutions of Type-I, and contain decelerating shock 
waves. New, Type-II, self-similar solutions for almost all the range u > 3 containing 
accelerating shock waves were constructed. 
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Here we briefly summarize the Type-II solutions for u > 3. The hydrodynamic 
equations for an ideal gas with adiabatic index 7 in spherical symmetry are given by: 

(d t + ud r )p + pr~ 2 d r (r 2 u) = , 
p(d t + ud T )u + d r (pc 2 /^) = , 



(1; 



(d t + ud r )(c 2 P l -y 1 ) = o , 

where the dependent variables u, c, and p are the fluid velocity, sound velocity, and density, 
respectively. We now seek a self-similar solution to the hydrodynamic equations (Eqn. [T]) of 
the form: 

u(r, t) = R£U(Z) , c(r,t) = R£C(0 , 

(2) 

p(r,t) = BR*G(0 , p(r,t) = BR*R 2 P(0 , 
where £ = r/R(t) is the dimensionless spatial coordinate, and the length scale R(t) 



(frequently abbreviated as simply R) is the shock radius and satisfies (Zel'Dovich and 
Raizer 1967 Waxman and Shvarts 1993[ ) 



RR 



6 



RocR 6 



(3) 



where 5 is a constant. The quantities G, C, U, and P, which are defined by Eqns. |2j give the 
spatial dependence of the hydrodynamic quantities. The diverging (exploding) solutions of 
Eqn. [3] are 

A(t-t ) a , 5<1 



S 



(4) 



A(t - t) a , 5 > 1 



where a = 1/(1 — 5). 



Solutions with 5 < 1 diverge in infinite time, and to represents the time of the point 
explosion, which is usually taken to be t = 0. For 5 < the shock wave decelerates and for 
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< 5 it accelerates. For 5 > 1 the shock wave accelerates so fast that it diverges in a finite 
time. In this case, to represents the time of divergence rather than the explosion time. The 
transition between finite and infinite divergence occurs at 5 = 1 where we have exponential 



time dependance (Simonsen and Meyer- Ter-Vehn 1997). 



Substituting Eqn. [2] into the hydrodynamic equations (Eqns. [T]) and using Eqn. |3j one 
gets regular differential equations for the similarity quantities U, C, and G (see for example 
Landau & Lifshitz) with two free constants, the similarity parameters e and 5: 

dU Ai (C/,C) dC A 2 {U,C) 



dlog£ A(U,C) ' dlog£ A(U,C) 
and an explicit expression for the density G: 



(5) 



C- 2 (l - u) x G~<- 1+x £, 3X ~ 2 = const . (6) 

The functions A, Ai, and A2 are given by: 

A = C 2 - (1 - Uf , 

Ai = U(l - U)(l -U -5)- 3UC 2 - 3C 2 (e + 25) h , 

(7) 

A 2 = C(l-U)(l-U-6)-(rf- l)CU{\ -U + 5/2)- 

3 1 25 - (7 - l)e 

and the parameter A is 

3 + e V ; 

The similarity parameter e can be found from the boundary conditions at the strong 



shock, the Hugoniot jump conditions (Landau and Lifshitz 1987). Applying these relations 



to a strong shock one gets e = —oj, and also 



urn = 4f , cm = , cm = 1+1 . ( 9) 

7 + 1 7 + 1 7-I 
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The boundary conditions on the shock do not state any limits on the possible values of 
the similarity parameter 5. In order to determine the value of this parameter one should 
distinguish two kinds of similarity flows: Type-I and Type-II, defined first by Zel'dovich 



(Zel'Dovich and Raizer 1967). A solution of Type-I describes the flow in all space and 



therefore conservation laws must be obeyed by the self-similar solution. One can then 
deduce 6 — (u — 3)/2, which gives the well-known Sedov- Taylor solutions. However, for 
u > 3 it is easy to see that the solution obtained with this value of 5 contains an infinite 
amount of energy and therefore can not describe the flow over the whole space. Therefore, 
the flow must be Type-II. 

In Type-II solutions, there is a region, whose scale relative to the flow characteristic 
length R(t) goes to zero with time, in which the similarity solution does not describe 
the physical system. Therefore, for this kind of solution the energy does not have to be 
conserved in the self-similar solution since this solution does not describe the whole flow. 
In order that the region which is not self-similar (located around the origin) does not 
influence the self-similar solution, the solution must pass through the singular point defined 



by (Zel'Dovich and Raizer 1967 Waxman and Shvarts 1993): 



U + C= 1 . 



(10) 



From this singular point requirement, the dependence of 5 upon the parameters u and 



7 can be found. It was found (Waxman and Shvarts 1993) that for oj > uj 9 (j) > 3 there is 
a value of 5 for which the solution passes through a singular point, and therefore a second 
type self-similar solution exists. 



A fully analytic solution to Eqns. [5]48] exists for the case where u> = co a (l) 
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2(47-l)/(7+l): 



7 + 1 ^ > ^ vs; 7 + l' 

(11) 

g(0 = ^jr 8 , p(o = ■ 

For this analytical case the parameter 5 is given by 5 = (7 — 1)7(7 + 1). 



3. Slightly Two- and Three-Dimensional Self-Similar Solutions 

We now consider small deviations from the spherically symmetric problem discussed 
in Sec. [2j In general, we wish to solve the problem for an external density perturbation of 
arbitrary angular dependence, which we parameterize by 



p{r,8,<f>) =r- u >[l + aF(9, 



(12) 



where a <C 1, F is an arbitrary function of 9 and 0, and uj > 3. 

However, before embarking on this general problem we assume in this and the following 
section that F(9,(f>) = Yi j7n . In Sec. [5] these solutions will be used to construct the solution 
for an arbitrary F(9, (ft). 

We shall use here the Eulerian perturbation approach. We define the perturbed 
quantities as the difference between the perturbed solution (i.e., the slightly two-dimensional 
self-similar solution) and the unperturbed one-dimensional solution at the same spatial 



point. The derivation of the perturbation equation is similar to the one given by Ryu and 



Vishniac (1987), Chevalier (1990), and Sari et al. (2000). The perturbed hydrodynamic 
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quantities are defined as 

Sv(r, 9, 0, t) = v{r, 9, 0, t) - v {r, t)f , 

Sp(r, 9, 0, t) = p(r, 9, 0, t) - p (r, t) , (13) 

5p(r, 9, 0, t) = p(r, 9, 0, t) - p (r, t) , 

where v, p, and p are the velocity, pressure, and density in the perturbed solution, while 
v r, p , and p are the same quantities as in the unperturbed solution. 



We consider perturbations that can be written in a separation of variables form (Cox 



1980): 



where 



Sv(r, 9, 0, t) = £R [8U r ^)Y lm {9, 0)f + 8U T {£)V T Y lm {6, 0)] / , 

8p(r, 9, 0, t) = BmG(0Y lm (9, 0)/ , (14) 
5p(r,9,(f>,t) = BR^R 2 6P(OY lm (0A)f , 



V T = 9^ + ^—J- (15) 
39 sinfcl ocb 



are the tangential components of the gradient and R(t) is the unperturbed shock radius 
which still satisfies Eqn. [3] The perturbed shock radius, R(t, 9, 0), is given by 



R(t, 9, 0) - R(t) = 6R(t, 9, 0) = Y lim (9, <f>)R{t)f . (16) 



Eqns. [141 and [TB] define the quantities 5U r , 5Ut, 8P, 5G, and /. The quantity / measures 
the fractional amplitude of the perturbation to the shock wave radius. Here we deviate 
from the standard treatment of stability. There, / is a function of time: if the function 
/ increases with time then the solution is unstable, while if / decreases with time then 
the solution is stable. However, here, since we demand that the perturbed solution be 
self-similar, / has to be independent of time. 
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We linearize the hydrodynamic equation around the unperturbed self-similar solution 
to get a linear set of equations: 

MY' = NY (17) 



where 



Y 



M 



5U T 



\ 6P J 



( 






G 



{u-\)ec 









(U-l)eG 




^ 

1 





N 



P'G^ 1 



V G 2 ' t{~P 
and G, U, and P are defined by Eqn. [2] 



-£G'-3G 1(1 + 1)G 

{l-6-2U-£U')G£ 
(l-6-U)G£ 










-r 1 
t(u-i; 



p2j 



Unlike the perturbation equations for stability, the equations above do not contain an 



unknown parameter. They are, in fact, a special case of the equations used in Sari et al. 



(2000), but with the perturbation growth rate set to q = 0. In that sense they are similar 



to the equations of Oren and Sari (2009) for discretely self-similar solutions. Instead, a 



new parameter d — a/f appears in the shock boundary conditions (the linearized Hugoniot 
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jump conditions): 

5G(l) = %±\(d-u)-G' , 6U r (l) = -U' 



7 



For any value of the parameter d one can integrate Eqn. 17 beginning at the shock 
front using the shock boundary conditions. However, the singular point of the unperturbed 
solution, £ c , where C + U = 1, is also a singular point of the perturbed solution. Therefore, 
in general, such integration will diverge at the sonic point £ c . Only for specific values of the 
parameter d, where an additional boundary condition at the singular point is satisfied, is 
the solution regular. These are the physical values for the parameter d. 

Technically, solving these equations is easier than the equivalent perturbation case. 
The reason is that the unknown parameter d appears only in the shock boundary condition, 
and is absent from the differential equations. We can therefore solve these equations 
starting from the sonic point outward, and find the three independent solutions that are 
nonsingular at £ c . Then we can find a linear combination of these three solutions, and the 
value of d that can solve the four boundary conditions at the shock. 



4. Results 

For convenience we investigate the case 7 = 5/3, u = 17/4, where the unperturbed 
solution is analytic. For I — 1 we obtain d = —11.2. This means that the fractional 
amplitude of perturbations in the shock wave position, /, are an order of magnitude smaller 
than the fractional amplitude of perturbations in the external density a . The negative 
sign implies that at angles where the external density is higher, the shock wave position is 
retarded. This is expected intuitively. From the shock boundary conditions, we infer that 
the pressure at these angles is also lower. For I = 2 we find d — —11.6, and for I = 3 we find 
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d — —12.1. A plot of d as function of I is given in Fig. [TJ 



5. Extension to Arbitrary Angular Dependence 



The analysis above is limited to external density perturbations whose angular 
dependence is a spherical harmonic. This is necessary in order to obtain separation 
between the angular and radial dependencies. However, since we are dealing with linear 
perturbations, any arbitrary angular dependence can be expanded into a sum of spherical 
harmonics, each of which could be solved in the method described in the previous section. 
Then, the solutions can be summed, leading to the perturbation solution for external 
density perturbations with arbitrary angular dependence. 

As an example, we consider the following problem: A strong point-like explosion is 
launched into a surrounding which has a density on one side of a plane slightly different 
from the density on the other side of the plane. In our notation this is p ex r~ 17 / 4 (l + oH{6)) 
where H(8) = 1 for 6 < n/2 and H{6) = —1 for 9 > tt/2. The point explosion in half space 
could be thought of as an extreme version of this density profile with a — 1. However, our 
solution formalism applies only for slightly two-dimensional cases where er <C 1. 

Such a density profile could be expanded in spherical harmonics as 




(19) 



The shape of the shock, R + 5R(8), deviates from its unperturbed value R by 




(20) 



This shape is plotted in Fig. [2] for both this analytic solution and for the numerical solution 
discussed in Sec. [7j To make the analytic curve the sum was taken from n = to n = 50. 



13 




20 40 60 80 



spherical harmonic degree 1 



Fig. 1. — Dots show d as function of I as obtained by solving the differential equations for 
7 = 5/3 and u = 17/4. For small I we have d = —10., while for large I, i.e., short wavelength, 
we obtain a linear relation: d = —\fhj\l (solid line). 
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Fig. 2. — The analytic (green) and numerical (blue) fractional deviation of the shock position 
as function of 9 for the Heaviside density distribution (Eqn. [l9| in units of a. The numerical 
result is obtained for o = 0.01, and is discussed in Sec. [7j while the analytic solution comes 



from Eqn. [20} It can be seen in both curves that, roughly speaking, the shock is composed 
of two hemispheres, connected smoothly over a short angular scale of less than 0.1 radians 
(FWHM). 
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6. Short Wavelength Limit 

Because the flow does not vary in the short wavelength limit, we may treat the matrices 
M and N as constants close to the shock front. By using the unperturbed values of the 
state variable at the shock, we find the four independent modes of the problem: 

A = -8,3,±y^L (21) 

The first two are independent of I and indicate that the state functions, close to the shock, 
vary on the scale R, regardless of the wavelength of the perturbation. However, the other 
two are linear in I meaning that close to the shock the state functions vary over small 
scales of order R/l. Therefore, for large /, the positive mode is growing inward very rapidly, 
and thus can not exist physically. For this reason we demand that the perturbation has 
no component along this mode on the shock front by requiring it to be written as a linear 
combination of the three eigenvectors associated with the other modes. This provides the 
extra boundary condition at the shock that allows us to determine d. Performing this 
calculation we find that for general u and 7 

d = -(wi l (22 » 

in the limit I ^> 1. We plot the general solution of d(l) in Fig. [I] for the case discussed in 
Sec. [4j along with the short wavelength limit described by Eqn. [22] It can be seen that the 
agreement is good. 



7. Comparison with 2D Numerical Hydrodynamical Simulations 



To compare the analytic solution presented in Sees. [2}j6] to numerical results we used 



the PLUTO hydrodynamic code (Mignone et al. 2007) to simulate an explosion on a weakly 



discontinuous surface with a = 0.01 (see Eqn. 12). Again for convenience we consider the 



case where the unperturbed solution is analytic: 00 = 17/4 and 7 = 5/3. 
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Fig. 3. — Comparison between the numerical (blue) and analytic (green) solutions for the 
normalized I = 1 self-similar pressure perturbation as a function of the self-similar variable £ 



for a Heaviside initial density distribution (Eqn. 19 ). Good agreement between the simulation 
and the analytic solution is found in the region of self-similarity (0.76 = £ c < £ < 1). Details 
of the simulation are discussed in Sec. [7j while the analytic solution is described in Sees. [2]-[6j 
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The computational mesh had 10 5 cells in the radial direction, and 100 cells in the 
tangential direction. The inner radius was 10 -6 and the outer radius was 1. The smallest 
angle was and the largest n. The radius of the initial hot spot was 2 x 10~ 4 , and the 
pressure there was 10 18 , whereas outside the hot spot the pressure was 1. The Riemann 
solver used was hllc. 



The numerical and analytic results are compared in Figs. [2}|6j Though the self-similar 
solution is valid everywhere, the deeper one looks into the flow, the longer it takes for 
the physical flow to approach this solution. Therefore, at any finite time, there exists an 
inner region that is not in agreement with the self-similar solution. Our comparison tends 
to reflect these points and in all cases there is agreement to within 10% throughout a 
significant fraction of the flow. 

In particular, in Fig. [2] we compare the deviation of the shock radius from the 
unperturbed solution to the shock radius in units of o in both cases. Analytically this 



function is independent of a (see Eqn. 20). Near the interface of the media (9 = vr/2) the 
two solutions match well, while near the poles (8 = 0, tt) the two differ by approximately 
10%. 

In the following four figures (Fig. |3j Fig. |4j Fig. [5j and Fig. [6]) we compare the fractional 
deviation of the pressure and angular velocity for both I = 1 and 1 = 2. As expected, as is 
the case with the general solution plotted in Fig. |2j discrepancies between the numerical 
and analytic work are always less than 10%. Descriptions of the specific cases are given in 
the captions. 
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Fig. 4. — Comparison between the numerical (blue) and analytic (green) solutions for the 
normalized I = 2 self-similar pressure perturbation as a function of the self-similar variable £ 



for a Heaviside initial density distribution (Eqn. 19 ). Good agreement between the simulation 
and the analytic solution is found in the region of self-similarity (0.76 = £ c < £ < 1). Details 
of the simulation are discussed in Sec. [7j while the analytic solution is described in Sees. [2]-[6j 



-19- 



X 10" 



I = 1 




-10 



0.5 0.55 0.6 



Fig. 5. — Comparison between the numerical (blue) and analytic (green) solutions for the 
normalized I = 1 self-similar fractional angular velocity as a function of the self-similar 



variable £ for a Heaviside initial density distribution (Eqn. 19). Good agreement between 



the simulation and the analytic solution is found even well outside the region of self-similarity 
(0.76 = £ c < £ < 1). Details of the simulation and its parameters are discussed in Sec. [7j 
while the analytic solution is described in Sees. pffSj 
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Fig. 6. — Comparison between the numerical (blue) and analytic (green) solutions for the 
normalized I = 2 self-similar fractional angular velocity as a function of the self-similar 



variable £ for a Heaviside initial density distribution (Eqn. 19). Good agreement between 
the simulation and the analytic solution is found even well outside the region of self-similarity 
(0.76 = £ c < £ < 1). Details of the simulation are discussed in Sec. [7| while the analytic 
solution is described in Sees. 



8. Discussion 



We have considered the problem of a strong shock propagating into a slightly aspherical 
medium made up of a density with a spherically symmetric radial power-law plus a 
perturbation of arbitrary angular dependence, and solved for the Type-II self-similar 
solution. Such an external medium has a density profile p(r, 6>, 0) = r _w [l + aF(Q,4>)\, 
where u > 3, c C 1, and F is an arbitrary function of 9 and 0. Because the perturbations 
are small, the hydrodynamic equations can be linearized around the unperturbed solution. 
This then allows us to expand F as a series in spherical harmonics, and solve the problem 
term by term. In this way the general problem is reduced to one which includes only 
perturbations F(6,<j)) oc Yi m (0,<j>). 

The linearized self-similar equations are presented for this simpler case, F(9, 0) oc 
Yi m (9,(p), along with the appropriate boundary conditions. There is a unique solution 
to these equations which depends on a single parameter d, which is determined by the 
requirement that the solution would pass smoothly through the sonic point. That the 
only dependence on d is in the boundary conditions makes these equations particularly 
straight-forward to solve. 

We demonstrate this process on a specific example which deviates from spherical 
symmetry by a weak step function in the outside density across a plane containing the 
initial explosion. As expected, instead of the shock being spherical, it is composed of 
two hemispheres smoothly connected across the plane of the discontinuity. Our 2D 
hydrodynamical simulations agree well with this solution. 

We thank Yonatan Oren for helpful discussions. This research was partially supported 
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